#ifndef MATRIX_ALGO
#define MATRIX_ALGO

#include "Matrix.h"
#include "Matrix_exception.h"
#include "Matrix_func_obj.h"

using namespace std;
using namespace Matrix_lib;


/** Multiply(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator*(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(Mul_assign<T>(),c);	return m; }

/** Add(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator+(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(Add_assign<T>(),c);	return m; }

/** Minus(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator-(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(Minus_assign<T>(),c);	return m; }

/** AND(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator&(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(And_assign<T>(),c);	return m; }

/** OR(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator|(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(Or_assign<T>(),c);	return m; }

/** XOR(Matrix<T,N>,scalar) operation **/
template<class T, int N>
Matrix<T,N> operator^(const Matrix<T,N>& m1, const T& c) { Matrix<T,N> m(m1); m.base_apply(Xor_assign<T>(),c);	return m; }


/** Apply function f on each element of calling matrix **/
template<class T, int N, class F> 
void apply(Matrix<T,N>& m, F f)            { m.base_apply(f); }

/** Apply function f with argument c on each element of calling matrix **/
template<class T, int N, class F> 
void apply(Matrix<T,N>& m, F f,const T& c) { m.base_apply(f,c); }
//arithmetic operations

/** Element-wise operation multiply operation on Matrix<T,N> elements **/
template<class T, int N>
void operator*=(Matrix<T,N>& m, const T& c) { m.base_apply(Mul_assign<T>(),c); }

/** Element-wise operation divide operation on Matrix<T,N> elements **/
template<class T, int N>
void operator/=(Matrix<T,N>& m, const T& c) { m.base_apply(Div_assign<T>(),c);  }

/** Element-wise operation modulus operation on Matrix<T,N> elements **/
template<class T, int N>
void operator%=(Matrix<T,N>& m, const T& c) { m.base_apply(Mod_assign<T>(),c); }

/** Element-wise operation multiply operation on Matrix<T,N> elements **/
template<class T, int N>
void operator+=(Matrix<T,N>& m, const T& c) { m.base_apply(Add_assign<T>(),c);  }

/** Element-wise operation minus operation on Matrix<T,N> elements **/
template<class T, int N>
void operator-=(Matrix<T,N>& m, const T& c) { m.base_apply(Minus_assign<T>(),c); }
//relational operations

/** Element-wise operation OR operation on Matrix<T,N> elements **/
template<class T, int N>
void operator|=(Matrix<T,N>& m, const T& c) { m.base_apply(Or_assign<T>(),c);   }

/** Element-wise operation XOR operation on Matrix<T,N> elements **/
template<class T, int N>
void operator^=(Matrix<T,N>& m, const T& c) { m.base_apply(Xor_assign<T>(),c); }

/** Element-wise operation AND operation on Matrix<T,N> elements **/
template<class T, int N>
void operator&=(Matrix<T,N>& m, const T& c) { m.base_apply(Xor_assign<T>(),c); }

/** Dot product of two vectors, return value is of type T **/
template<class T>
T& sdot_mat(Matrix<T,1>& v1, Matrix<T,1>& v2) {
	T result = 0;
	if(!v1.size() || !v2.size() || (v1.dim1() != v2.dim1())) throw Matrix_exception("Invalid matrix sizes for sdot_mat");
	Index dim = v1.dim1();
	for(Index i=0; i<dim; ++i)
		result += v1(i)*v2(i);
	return result;
}

/** Swapping two Matrices, both matrices should be Matrix<T,N> **/
template<class T, int N>
void swap_mat(Matrix<T,N>& m1, Matrix<T,N>& m2) {
	if(!m1.size() || !m2.size()) throw Matrix_exception("Invalid matrix sizes for swap_mat");
	Matrix<T,N> temp = m1;
	m1=m2;
	m2=temp;
}

/** Operation is of form y = y + alpha*x, where x & y are Matrices of same size and dimensions **/
template<class T, int N>
Matrix<T,N> scale_and_copy(const Matrix<T,N>& m1, const Matrix<T,N>& m2, const T& x) {
	if(!m1.size() || !m2.size()) throw Matrix_exception("Invalid matrix sizes for scale_and_copy");
	return (m1 + m2*x);
}

/** Matrix<T,2> * Matrix<T,1> product with Matrix<T,1> as result **/
template<class T>
Matrix<T,1>& prod_mv(Matrix<T,2>& m, const Matrix<T,1>& v) {
	if(!m.size() || !v.size()) throw Matrix_exception("Invalid matrix sizes");
	if(m.dim2() != v.dim1()) throw Matrix_exception("Matrix dimensions don't match for prod_mv operation");
	Matrix<T,1> temp(m.dim1());
	for(Index i=0;i<m.dim1();++i)
		for(Index j=0;j<m.dim2();++j)
			temp(i)=m(i,j)*v(j);
	return temp;
}

/** Transpose of Matrix<T,1>  **/
template<class T>
Matrix<T,2>& vector_transpose(const Matrix<T,1>& v) {
	if(!v.size()) throw Matrix_exception("Invalid vector size vector_transpose");
	return (Matrix<T,2>(v));
}

/** Transpose of Matrix<T,2>  **/
template<class T>
Matrix<T,2>& matrix_transpose(Matrix<T,2>& m) {
	if(!m.size()) throw Matrix_exception("Invalid matrix size matrix_transpose");
	Matrix<T,2> temp(m.dim1(), m.dim2());
	for(Index i=0;i<m.dim1();++i)
		for(Index j=0;j<m.dim2();++j)
			temp(i,j)=m(j,i);
	return temp;
}

// /**			Level 3			**/
//Level 3 operations can be defined using MAtrix class member operations and operations given above


#endif

